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Report  Title 

Homogeneous  Self-Dual  Algorithms  for  Stochastic  Semidefmite  Programming 

ABSTRACT 

Ariyawansa  and  Zhu  [3]  have  proposed  a  new  class  of  optimization  problems  termed 

stochastic  semidefmite  programs  (SSDPs)  to  handle  data  uncertainty  in  applications  leading  to  (deterministic) 
semidefmite  programs  (DSDPs).  For  the  case  where  the  event  space  of  the  random  variables  in  an  SSDP  is  finite  they 
have  also  derived  a  classof  volumetric  barrier  decomposition  algorithms,  and  proved  polynomial  complexityof  the 
short-step  and  long-step  members  of  the  class  [2],  When  the  event  space  of  therandom  variables  in  an  SSDP  is  finite, 
the  SSDP  is  equivalent  to  a  large  scale  DSDPwith  special  structure.  Polynomial  homogeneous  self-dual  algorithms  [11] 
are  animportant  class  of  algorithms  that  have  been  introduced  for  solving  (general)  DSDPs.lt  is  therefore  possible  to 
solve  SSDPs  by  applying  homogeneous  self-dual  algorithmsto  their  DSDP  equivalents.  However,  such  algorithms, 
while  polynomial,  will  still 

have  high  computational  complexities  in  comparison  to  decomposition  algorithms. 

In  this  paper,  we  show  how  the  special  structure  in  DSDP  equivalents  of  SSDPs 
can  be  exploited  to  design  homogeneous  self-dual  algorithms  with  computational 
complexities  similar  to  those  of  volumetric  barrier  decomposition  algorithms. 


Homogeneous  Self-Dual  Algorithms  for  Stochastic 
Semidehnite  Programming 

Shengping  Jin*  •  K.  A.  Ariyawansa^  •  Yuntao  Zhih 


Abstract 

Ariyawansa  and  Zhu  [3]  have  proposed  a  new  class  of  optimization  problems  termed 
stochastic  semidehnite  programs  (SSDPs)  to  handle  data  uncertainty  in  applications 
leading  to  (deterministic)  semidehnite  programs  (DSDPs).  For  the  case  where  the 
event  space  of  the  random  variables  in  an  SSDP  is  finite  they  have  also  derived  a  class 
of  volumetric  barrier  decomposition  algorithms,  and  proved  polynomial  complexity 
of  the  short-step  and  long-step  members  of  the  class  [2] .  When  the  event  space  of  the 
random  variables  in  an  SSDP  is  finite,  the  SSDP  is  equivalent  to  a  large  scale  DSDP 
with  special  structure.  Polynomial  homogeneous  self-dual  algorithms  [11]  are  an 
important  class  of  algorithms  that  have  been  introduced  for  solving  (general)  DSDPs. 
It  is  therefore  possible  to  solve  SSDPs  by  applying  homogeneous  self-dual  algorithms 
to  their  DSDP  equivalents.  However,  such  algorithms,  while  polynomial,  will  still 
have  high  computational  complexities  in  comparison  to  decomposition  algorithms. 
In  this  paper,  we  show  how  the  special  structure  in  DSDP  equivalents  of  SSDPs 
can  be  exploited  to  design  homogeneous  self-dual  algorithms  with  computational 
complexities  similar  to  those  of  volumetric  barrier  decomposition  algorithms. 

Keywords:  Semidehnite  programming,  homogeneous  self-dual  algorithms,  compu¬ 
tational  complexity,  stochastic  semidehnite  programming. 


1  Introduction 

Stochastic  programming  models  have  been  very  useful  in  dealing  with  data  uncertainty 
in  many  applications.  Stochastic  linear  programs  were  introduced  in  the  1950s  as  a 
paradigm  for  dealing  with  uncertainty  associated  with  data  in  applications  leading  to 
linear  programs.  Since  then  they  have  been  studied  extensively  [6,  8,  18,  19]. 

Deterministic  semidehnite  programs  (DSDPs)  have  been  the  focus  of  intense  research 
during  the  past  20  years,  especially  in  the  context  of  interior  point  methods  for  opti- 
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mization  [1,  10,  13,  16].  DSDPs  generalize  linear  programs  and  have  a  wide  variety  of 
applications,  especially  beyond  those  covered  by  linear  programs. 

More  recently,  Ariyawansa  and  Zhu  [3]  (see  also  [9])  introduced  a  new  model  termed 
Stochastic  Semidefinite  Programs  (SSDPs)  to  handle  data  uncertainty  in  applications 
leading  to  DSDPs.  In  this  new  model,  they  combine  stochastic  programming  and 
semidefinite  programming  together  so  that  the  uncertainty  in  data  leading  to  a  semidefi¬ 
nite  program  can  be  dealt  with  in  the  same  way  that  a  stochastic  linear  program  handles 
uncertainty  in  data  in  an  application  leading  to  a  linear  program.  At  the  same  time,  since 
semidefinite  programs  generalize  linear  programs,  stochastic  semidefinite  programs  are 
also  a  generalization  of  stochastic  linear  programs.  See  [21]  for  a  preliminary  collection 
of  applications  of  SSDPs. 

Ariyawansa  and  Zhu  [2]  (see  also  [9])  have  presented  a  class  of  volumetric  barrier 
decomposition  algorithms  for  SSDPs  and  proved  polynomial  complexity  of  the  short- 
step  and  long-step  members  of  the  class.  Their  derivation  is  for  the  case  where  the  event 
space  of  the  random  variables  in  the  SSDP  is  finite.  The  computational  complexity  in 
terms  of  the  number  of  arithmetic  operations  of  the  short-step  and  long-step  algorithms 
in  [2]  are  0(K 1'5)  and  0(K2)  respectively,  where  K  is  the  number  of  realizations  of  the 
random  variables.  Another  important  feature  of  the  algorithms  in  [2]  is  that  the  most 
expensive  part  of  the  algorithm  naturally  separates  into  K  subproblems,  which  allows 
efficient  parallel  implementations. 

When  the  event  space  of  the  random  variables  in  an  SSDP  is  finite,  the  SSDP  can 
be  converted  to  an  equivalent  large  scale  DSDP.  Thus  such  an  SSDP  can  also  be  solved 
by  using  algorithms  for  DSDPs  on  the  DSDP  equivalent  to  the  SSDP.  One  such  general 
purpose  DSDP  algorithms  is  the  homogeneous  self-dual  algorithm  (see  [11]).  However, 
the  computational  complexity  in  terms  of  the  number  of  arithmetic  operations  of  the 
general-purpose  homogeneous  self-dual  algorithm  on  the  DSDP  equivalent  of  an  SSDP 
is  0(K4  5).  In  the  context  of  SSDPs,  K  is  large,  and  so  0(/P4"5)  is  significantly  larger 
than  0(K 1"5)  and  0(K 2),  the  complexities  of  the  short-step  and  long-step  algorithms 
respectively  in  [2], 

In  this  paper,  we  propose  homogeneous  self-dual  algorithms  for  SSDPs  (with  finite 
event  spaces  of  their  random  variables)  that  exploit  the  special  structure  of  the  equiva¬ 
lent  DSDP  to  reduce  the  computational  complexity  to  0(K 1"5).  In  addition,  the  most 
expensive  part  of  the  algorithm  separates  into  K  subproblems. 

In  general,  algorithms  for  stochastic  optimization  belong  to  two  broad  categories  both 
exploiting  their  special  structure  but  in  two  different  ways.  First,  we  have  algorithms 
based  on  notions  of  cutting  planes  and  decomposition.  The  classical  L-shaped  algorithm 
in  [17],  the  volumetric  center  cutting  plane  algorithm  in  [4],  and  the  decomposition 
algorithms  in  [20],  [9]  and  [2]  belong  to  this  category.  Second,  we  have  algorithms 
designed  by  tailoring  the  operations  in  general  purpose  algorithms  to  exploit  the  special 
structure  in  stochastic  optimization  problems.  The  algorithms  in  [12,  7,  5]  belong  to 
this  category.  The  algorithms  derived  in  this  paper  for  SSDPs  also  belong  to  this  second 
category. 

The  rest  of  the  paper  is  organized  as  follows.  In  the  next  section  we  introduce  our 
notation  and  present  some  preliminary  material  on  the  homogeneous  self-dual  methods 
for  solving  DSDPs.  We  present  our  homogeneous  self-dual  algorithms  for  SSDPs  in  §3. 
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In  §4  we  present  an  efficient  method  that  exploits  the  structure  of  the  DSDP  equivalent 
of  an  SSDP  for  computing  search  directions  in  homogeneous  self-dual  methods.  We  also 
obtain  an  estimate  of  the  number  of  arithmetic  operations  per  iteration  of  our  algorithm 
in  §4.  We  conclude  the  paper  in  §5  where  we  compare  the  computational  complexities 
of  the  new  algorithms  presented  in  this  paper  with  those  of  the  algorithms  in  [2], 


2  Preliminaries 


2.1  Notation  and  terminology 

The  notation  and  terminology  we  use  in  the  rest  of  this  paper  follow  that  in  [13].  We 
let  Rmxn  and  MnVn  denote  the  vector  spaces  of  real  mxn  matrices  and  real  symmetric 
n  x  n  matrices  respectively.  We  write  X  P  0  ( X  >-  0)  to  mean  that  X  is  positive 
semidefinite  (positive  definite)  and  we  use  U  P  V  or  V  A  U  to  mean  that  U  —  V  P  0. 
For  U,  V  £  Mmxn  we  write  U  •  V  :=  trace  (UT V)  to  denote  the  Frobenius  inner  product 
of  U  and  V. 

Given  At  £  MnVn  for  i  =  1,  2, . . . ,  m,  we  define  the  operator  A  :  MnVn  — >  Mm  by 


AX 


Ai»X 

a2»x 

Am  •  X 


(1) 


for  any  X  £  MnVn. 

The  adjoint  operator  of  A  with  respect  to  the  standard  inner  products  in  MnVn  and 
Mm  is  the  operator  A*  :  Km  — >  MnVn  defined  by 


A*y  =  ^2viAi,  (2) 

i= 1 


for  any  y  £  Mm. 

We  introduce  the  useful  notation  Kronecker  product  P  ®  Q,  where  P  and  Q  are 
usually  symmetric.  This  is  an  operator  from  MnVn  to  itself  defined  by 

(P®Q)U  =^(PUQJ +  QUPJ),  (3) 

for  any  U  £  MnVn. 

Finally  we  introduce  the  operator  svec  which  transforms  a  matrix  in  MnVn  to  a  vector 
in  M”  (see  [14])  where  n  :=  \n{n  +  1). 


2.2  Homogeneous  Self-Dual  Methods  for  DSDPs 

Given  data  d*  £  MnVn  for  i  =  1,  2, . . . ,  m,  b  £  Rm  and  C  £  RnVn,  a  DSDP  in  primal 
standard  form  is  defined  as 

minimize  C  •  X 

subject  to  AX  =  b  (4) 

X  h  0, 
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where  X  G  MnVn  is  the  variable.  The  dual  of  (4)  is 

maximize  bTy 

subject  to  A*y  +  S  =  C  (5) 

S  ^  0, 


where  y  G  and  S  G  MnVn  are  the  variables. 

We  briefly  review  the  homogeneous  interior  point  algorithm  for  DSDPs  as  in  [11], 
The  homogeneous  model  for  (4,  5)  is  as  follows: 

AX  -br  =0 

-A*y  -S  +tC  =  0  (6) 

-C  •  X  +bJy  — k  =  0. 

It  is  easy  to  show  from  (6)  that 

X  •  S  +  tk  =  0. 

The  main  step  at  each  iteration  of  the  homogeneous  interior  point  algorithm  (shown 
below  in  Algorithm  1)  is  the  computation  of  the  search  direction  (AX,Ay,AZ)  from 
the  symmetrized  Newton  equations  with  respect  to  an  invertible  matrix  P  (which  is 
chosen  as  a  function  of  ( X ,  y,  S ))  defined  by: 


AAX 

—bAr 

=  rjrp 

-A*  Ay 

-AS 

+AtC 

=  vRd 

•  AX 

+bJ  Ay 

—A  k  =  rjrg 

kAt 

+tAk  =  7  n  —  TK 

£AX 

+PAS 

=  7  ^I-HP{XS) 

where  rp  :=  br—AX,  Rd  :=  A*y+S— tC,  rg  :=  C*X— bJy+K,  fi  :=  [1/(77.+  1)](A»5+tk), 
y  and  7  are  two  parameters,  Hp  :  MnVn  — >-  MnVn  is  the  symmetrization  operator  defined 
by 

Hp{U)  :=  l~{PUP~l  +  (P~1)TUJPT), 

and  £  :  MnVn  — ■>  MnVn  and  T  :  MnVn  — ■>  MnVn  are  the  linear  operators  defined  by 

£  :=P©(P"1)t5;  T  :=  PX  ©  (P”1)7, 


respectively. 

Currently,  the  most  common  choices  of  symmetrization  for  search  directions  in  prac¬ 
tice  are  as  follows  [15]. 

1.  Helmberg-Rendel-Vanderbei-Wolkowicz/Kojima-Shindoh-Hara/Monteiro  (HKM)  di¬ 
rection,  corresponding  to  P  :=  S'1/2.  In  this  case,  we  have  that  £  =  X  and 
P  =  X®S~1. 

2.  Kojima,  Shindoh,  and  Hara  (KSH)  direction  (rediscovered  by  Monteiro),  corre¬ 
sponding  to  P  :=  A-1/2.  In  this  case,  we  have  that  £  =  S  ©  X^1  and  T  =  X. 
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3.  Nesterov-Todd  (NT)  direction,  corresponding  to  P  :=  H  1//2,  here  H  is  the  unique 
symmetric  positive  definite  matrix  satisfying  HXH  =  X,  which  can  be  calculated 
by 

H  =  X1/2  (x1/2SX1/2y1/2  X1/2. 

In  this  case,  we  have  £  =  X  and  T  =  H  ©  H . 

Lemma  1.  Suppose  that  X  y  0,  S  >-  0,  and  A\,  A2, . . . ,  Am  are  linearly  independent. 
Then  for  each  of  the  above  three  choices  of  P,  £_1  and  exist,  and  £~XT  and  P~l£ 
are  positive  definite  and  self-adjoint. 

Proof.  See  [13].  □ 

We  state  the  generic  homogeneous  algorithm  as  in  [11], 


Algorithm  1  Generic  Homogeneous  Self-Dual  Algorithm  for  Solving  (4,  5) 
( X,y,S,T,K )  :=  (1,0,7, 1, 1) 
while  a  stopping  criterion  is  not  satisfied  do 
choose  i) ,  7 

compute  the  solution  (AT,  Ay,  AS,  At,  Ak)  of  the  linear  system  (7) 
compute  a  step  length  9  so  that 
X  +  6 AX  >-  0, 

S  +  6AS  P  0, 
r  +  6 At  >  0,  and 
k  +  9 An  >  0 

(A,  y ,  S,  t,  k)  :=  (X,  y,  S,  r,  k)  +  9 (AX,  Ay,  AS,  At,  Ak) 

end  while 


3  Homogeneous  Self-Dual  Algorithms  for  SSDPs 

3.1  Definition  of  a  SSDP  in  Primal  Standard  Form 

As  in  [3],  we  define  a  SSDP  with  recourse  in  primal  standard  form  based  on  deterministic 
data  Woi  G  ]RnoVno  for  i  =  1,2, . . . ,  mo,  ho  G  Mm°  and  Co  G  l,loVn0;  and  random  data 
Bi(u)  GK"»Vn»,  Wi(u)  G  MniVni  for  i  =  1,2,..., mi,  h(u)  G  Mmi  and  C(u)  G  M”iVni 
that  depend  on  an  underlying  outcome  u  in  an  event  space  P  with  a  known  probability 
function  P.  Given  this  data,  a  SSDP  with  recourse  in  primal  standard  form  is 

minimize  Co  •  Xq  +  E  [Q  (Xo,  w)] 

subject  to  yV(]X0  =  ho  (8) 

Xq  t  0, 

where  Xq  G  Mr!°Vno  is  the  first-stage  decision  variable,  Q  (Aq,  oj)  is  the  minimum  of  the 
problem 

minimize  C(uj)  •  Y (10) 

subject  to  B(co)Xo  +  W(w)Y(u)  =  h(u)  (9) 

Y(J)  h  0, 
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where  Y (uj)  G  M”lVni  is  the  second-stage  variable,  and 

E[Q(X0,u)]=  [  Q(X0,u)F(dw).  (10) 

Jn 

When  the  event  space  Q  is  finite,  Problem  (8,  9,  10)  is  equivalent  to  a  DSDP.  The 
equivalent  DSDP  has  a  special  block  angular  structure  which  is  exploited  to  reduce  the 
computational  complexity. 


3.2  Solving  SSDPs  with  Finite  Event  Space 

We  consider  SSDP  (8,  9,  10)  when  the  event  space  P  is  finite.  Let  {(Ck,  £y;.  V\4,  %),  k  = 
1,2,  ...,K}  be  the  possible  values  of  the  random  variables  (C(u),B(ui),W(uj),h(u)) 
and  let  pk  be  the  corresponding  probabilities  for  k  =  1,2, .. . ,  K.  For  convenience,  let 
Ck  '■=  PkCk  for  k  =  1, 2, . . . ,  K .  Then  Problem  (8,  9,  10)  is  equivalent  to 


minimize  Co  •  Xq  + 

C  i.Xi  +  • 

•  •  +  Ck  •  XK 

subject  to  Wo^fo 

=  ho 

BiX0  + 

WiXi 

=  hi 

BkX0 

+  WrXk  =  hx 

X0, 

Vi, 

o 

AI 

where  Xq  G  Mn°Vn°  and  Xk  G  M”lVni  for  k  =  1,  2, . . . ,  K  are  the  first-stage  and  second- 
stage  variables  respectively. 

Problem  (11)  is  a  DSDP  in  primal  standard  form.  To  see  this,  first  let 


C  :=  diag(C0,  Cu...,  CK)  G  R^o+^mMno+Kn^ 
X  :=  diag(X 0,  Xu...,  XK)  G  R("o+ifni)v(n0+Kni); 


and 


hJ  :  = 


Next  for  i  =  1,  2, . . . ,  (mo  +  Kmi),  dehne  Ai  G  M.(no+Kni)v(no+Kni^  by 


Ai  :=  diag(WOi,0, ...  ,0),  %  =  1,  2, . . .  ,m0 

Amo+i  ■=  diag{B\i,  Wi{,0, ...  ,0),  i  =  1,2, . . .  ,mi 

i)mi+i  ■  diag^Bxi  ,0, . . .  ,0, 1 Tjq  ),  i  1,2,...,  tti\  . 

Then  (11)  becomes 

minimize  C  •  X 

subject  to  Ai  •  X  =  hi,  i  =  1, 2, . . . ,  (mo  +  Km{)  (12) 

X  y  0. 


Problem  (12)  is  same  as  (4)  if  we  set  n  :=  (no  +  Kn\)  and  m  :=  (mo  +  Km\). 
So  Problem  (11)  is  a  DSDP  in  primal  standard  form  with  block  diagonal  structure. 
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Algorithms  that  exploit  this  special  structure  is  especially  important  when  K  is  large  as 
is  the  case  in  typical  applications. 

The  dual  of  (11)  is 


maximize  /ij  yo 

+  hjyi  + 

■  ■  +  hTKyK 

subject  to  W*yo 

+  B\yv  +  • 

■  ■  +  B*KyK 

+ 

So  = 

Co 

W\yx 

+ 

Si  = 

Ci 

W*KyK 

+ 

Sk  = 

Ck 

So , 

Si, 

Sr 

y 

o, 

where  y  G  R(,no+Ji'™i)  and  sk  g  flj(no+Ani)v(no+Ani)  for  k  =  1,2, . . . ,  K  are  the  variables. 

Now  the  homogeneous  interior  point  method  can  be  applied  to  Problem  (11,  13). 
However,  the  size  of  Problem  (11,  13)  increases  as  K  increases.  In  practice  K  is  typically 
very  large.  So  in  practice,  the  Problem  (11,  13)  is  large-scale  and  the  computation  of  the 
search  direction  in  Algorithm  1  (i.e.  the  solution  of  the  system  7)  is  very  expensive.  As 
we  shall  see  in  the  next  section,  this  computational  work  can  be  reduced  significantly  by 
exploiting  the  special  structure  of  Problem  (11,  13).  The  method  we  describe  in  the  next 
section  for  the  computation  of  the  search  direction  has  an  additional  desirable  feature: 
it  decomposes  into  K  smaller  computations  that  can  be  performed  in  parallel. 


4  An  Efficient  Method  for  Computing  Search  Directions 


We  now  describe  a  method  for  computing  the  search  direction  in  Algorithm  1  that 
exploits  the  special  structure  in  (11,  13). 

The  homogeneous  model  (6)  for  Problem(ll,  13)  is 


W0Ao  -  h0T 
BkXo  +  WkXk  —  hkT 

-  W0*y0  -  Ef=  i  &hVk  +  rC0  -  So 

-W*kyk  +  rCk  -  Sk 

Efc=0  hjyk  ~  E/c= o  Cc  •  xk  —  k 

Xk  >:  0,  Sk  >z  0,  k  =  0, 1,  2, ... ,  K,  r  >  0,  k  >  0, 


0 

0,  k  =  1,2,...,  K 

°  (14) 
0,  k  =  1,2,..., a:  1  f 

0 


where  Wk  and  B*k  for  k  =  1,2 are  adjoint  operators  in  the  sense  of  (2)  with 
appropriate  dimensions. 

The  search  direction  system  corresponding  to  (14)  can  be  derived  via  (7)  as 


W0AA0  -  Ji0At 
BkAXo  +  WkAXk  —  hkAr 

— Wq  Ay0  —  Efc=i  BkAyk+ 

AtCo  -  ASo 
-W*kAyk  +  A rCk  -  A Sk 
£0AXo  +  T0ASo 
£kXXk  +  XkA  Sk 
kAt  +  tAk 

Efc= o  hlAyk  -  Efc=o  ck  •  Axk  —  Ak 


yrpo 

r]TPk,  k  =  l,2,...,K 

vRdo 

T]Rdk,  k  =  1,2,. ..  ,K  (15) 

Ikh  ~  HPo(X0So) 

7 kh.  ~  HPk  (XkSk),  k  =  1, 2, . . . ,  K 
-yy  —  tk 

Vrg , 
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where 


rp o  =  h0T  -  W0X0 

rPk  =  hkr  -  BkX o  -  WfcXfe,  A;  =  1,  2, . . . ,  K 

Rd0  =  WqUo  +  EfcLl  BkDk  +  Sq  —  tCo 

Rdk  =  w kUk  +  sk  -  rCk,  k  =  1,2,...,  K 

r9  =  K  -  Ef=0  K'Vk  +  ELo  °k  •  Xk 

B  =  no+Km+1  (Ek= 0  *0>So  +  Tk), 

7]  and  7  are  two  parameters,  and  £k,  J~k  and  HPk  are  linear  operators  which  depend  only 
on  Xk  and  Sk . 

Now  we  present  the  crux  of  our  method  for  finding  the  search  direction  as  a  solution 
to  (15).  By  the  sixth  equation  of  (15)  and  the  fact  that  for  k  =  1,2, ...  ,K  exist, 
we  have 


A Sk  =  -Tkx£kA Xk  +  J^\ltiIk-HPk(XkSk)),  k  =  1, 2, . . . ,  K.  (16) 

Substituting  (16)  into  the  fourth  equation  of  (15),  we  get 

-W*kAyk  +  ArCk  +  Xkl£kAXk  =  r,Rdk  +  (7  -  HPk  (. XkSk ))  (17) 

for  k  =  1, 2, . . . ,  K .  Since  £j71  exist,  we  have  (Jrj71£k)1  =  £k  1  3~k ,  for  k  =  1,  2, . . . ,  K. 
By  equation  (17),  for  k  =  1,  2, ... ,  K ,  AXk  can  be  expressed  as 

=  S^FkWZAyk  -  S^MArCk  +  r,Rdk)  +  S^nh  -  HPk (. XkSk )).  (18) 

Substituting  (18)  into  the  second  equation  of  (15),  we  have 

Bk  AX0  +  m^XkWlAyk  -  £-lFk{ArCk  +  VRdk )  +  £k\^Ik  -  HPk(XkSk ))) 
-hkAr  =  rirPk,  k  =  1,2, . . .  ,K. 

(19) 

Thus 

Ayfc  =  MklBk  AX0  +  qkAr  +  vk  (20) 


where 


Mk  =  Wk£^rkW*k 

qk  =  M-x{Wk£-xTk  +  hk)  (21) 

uk  =  M~l(r]rPk  -  ■qWk£k1FkRdk  -  Wk£^l{^ylk  -  HPk(XkSk ))), 


for  k  =  1, 2, . . . ,  K . 

By  the  fifth  equation  of  (15)  and  the  fact  that  Rq1  exists,  one  gets 

A  So  =  -^SqAXq  +  -  HPo(X0S0)).  (22) 


We  substitute  (20)  and  (22)  in  the  third  equation  of  (15)  to  get 

— W0*Ay0  -  Ef=i  B%{-M~lBk AX0  +  qk At  +  uk)  +  AtCq  +  ^SqAJCq 
-  HPo(X0S0))  =  VRdo. 


(23) 


From  (23)  we  have 


AX0  =  ^o1(W0*Aj,0  +  (Ef=i^gfc-C'o)Ar) 

+M0" l(Ef=i  ^  +  vRd0  +  40  Wo  -  HPo(X0S0 )))  (24) 

=  M0- ‘W0*  Ayo  -  TqAt  +  UQ, 

where 

A4o  =  ^o  +  ELW1^- 

T0  =  WWEtiW  (25) 

4)  =  -Mo_1(Ef=i  +  r/^o  +  40  Wo  -  HPo(X0S0 ))). 

Now  substituting  (24)  into  the  first  equation  of  (15),  we  have 

W0(A4q^ lW(jAy0  -  T0  +  U0)  -  Ji0At  =  r]rPo.  (26) 


From  (26)  and  the  fact  that  WoAi0  1W g  is  nonsingular  (this  will  be  discussed  in  detail 
later  in  §5),  we  have  that 


Ay0  =  (WoA4o1Wo)_1((WoTo  +  /io)Ar  +  7?rpo  +  WoC/o 

=  a0Ar  +  /30, 


where 

«o  =  (WoM^WSyHWoTo  +  ho) 

Po  =  (W0W>W Wo  ~  ^Wo)- 

Now  we  substitute  backwards.  First  we  substitute  (27)  in  (24)  to  get 

AX0  =  Mo1W0*(aoAr  +  /3o)-ToAr  +  [/o 
=  f0Ar  +  $0, 

where 

*o  =  M^WSoq-Tq 
=  Mq'WSPq  +  Uq. 

Substituting  (29)  in  (20),  we  obtain 

A?/i  =  -(Wfcf’Ar1Tfc>Vfc)_113fc(4'oAr  +  4>0)  +  qkAr  + 
=  afcAr  +  (3k, 


(27) 

(28) 

(29) 

(30) 

(31) 


where 

afc  =  -(V^WW'Wo  +  gfc 

Pk  =  -(Wk£k1:FkW*k))-1)3k®o  +  i'k 

for  k  =  1, 2, . . . ,  X.  Furthermore,  we  substitute  (31)  in  (18)  to  get 


(32) 


AX, 


where 


^^fcWfcKAr  +  AO  -  £pl  Xk(ArCk  +  r?i74)  +  £^(7^/*  -  HPk(XkSk)) 

^fcA  r  +  ^fc,  k  =  1,2,...,X, 

(33) 

^4  =  4  —  A..  1RkCk  . 

<4  =  T01X,W,!/3fe-T01^?^+^71(7T4-^pfc(Afc5fc)).  1  J 
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for  A:  =  1,2,...,  K.  Now,  we  substitute  (27),  (29),  (31)  and  (33)  in  the  last  equation  of 
(15).  By  the  seventh  equation  of  (15)  this  yields 

K  K 

y  hJ(akAT  +  Pk)  -  y  Ck  •  (’LfcAr  +  <t>fc) - {-kAt  +  7 n~  tk)  =  r/rg. 

k= 0  k=0 


Finally  At  is  given  by 

Ar  =  TWg  +  r  J2k=0(Ck  •Qk-  h[pk)  +  (7 h  ~  tk)  5 

T  Ylk=o(hkak  ~  Ck  •  ^k)  +  « 

All  the  other  directions  can  be  obtained  by  (35).  Once  all  the  directions  are  computed, 
we  can  iterate  to  the  next  point  in  Algorithm  1. 


5  Complexity  Analysis 

In  this  section  we  first  show  that  under  reasonable  conditions  the  operations  described 
above  are  valid.  Then  we  estimate  the  computational  complexity  of  Algorithm  1  with  the 
method  described  in  §4  applied  on  Problem  (11,  13).  Finally,  we  compare  and  compare 
that  complexity  to  the  complexity  of  Algorithm  1  applied  on  Problem  (11,  13)  treating 
it  as  a  generic  primal-dual  DSDP  pair. 

5.1  Validation  of  Computations 

We  assume  that  IT0(1),  W0(2), . . . ,  TC0(mo),  and  wf w£2) . . . ,  W^mi)  for  k  =  1,  2, . . . ,  K  are 
linearly  independent.  Then  Mk  for  k  =  1,2 , ,K  in  21)  are  nonsingular  and  positive 
definite  by  Lemma  1. 

Now  we  will  show  that  Ato  in  (25)  is  nonsingular  and  that  tVo-Adg  1Wg  is  also  non¬ 
singular. 

Lemma  2.  Suppose  that  W^l\  Wq2\  . . . ,  IFq”!°^  and  that  W^\  . . . ,  Wk'1^  are  lin¬ 

early  independent  for  k  =  1,2, ,  K .  Then  A4q  in  (25)  and  Wo-Mg  1W’g  are  positive 
definite. 

Proof.  From  (25)  we  have 

Mo  =  ^-%  +  Ef.i  WB*. 

We  have  that  Jr(f1£o  is  positive  definite  by  Lemma  1,  so  it  suffices  to  show  that  BkMkxBk 
is  positive  semidehnite  for  k  =  1,2, ... ,  K .  In  fact,  denoting 

BkU  =  [b^  •  U,  B{2)  •[/,...,  B £mi)  •  u\ T 

for  each  U  G  M”°Vn°  and  Mjf1  =  ,  we  have 

L  J  J  m\  xmi 

B*kM^BkU  =  B*kW  1  \b^.U,...,B^.uY 

L  J  JmiXmi  L  J 
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for  k  =  1,  2, ... ,  K .  So 

BlM^BkU.U  =  T.iT.M''B'i)'U)l t'U 

=  (BkU)T  M~\BkU)  >  0. 

The  last  inequality  is  due  to  the  fact  that  Mr1  is  positive  definite.  As  in  [13]  we  can 
then  show  that  Wo-Mq  1Wg  is  positive  definite.  This  completes  the  proof.  □ 

5.2  The  Complexity  Analysis 

Theorem  1.  Suppose  that  m\  =  mo,  n\  =  no  and  that  wj^\  wj.2\  . . . ,  are  lin¬ 

early  independent  for  k  =  0, 1, . . . ,  K.  By  utilizing  method  described  in  for  computing 
the  search  directions  in  Algorithm  1,  we  have  that  the  number  of  arithmetic  operations 
in  each  iteration  of  Algorithm  1  is  0(K(m ,jj  +  m^n^)  +  n(j). 

Proof.  The  dominant  computations  of  the  method  in  §4  occur  at  (21),  (25),  (28),  (30), 
(32)  and  (34).  The  corresponding  numbers  of  arithmetic  operations  of  these  computa¬ 
tions  are  listed  in  Table  1.  In  particular,  the  computation  in  (25)  will  be  analyzed  in  de¬ 
tail.  The  total  number  of  arithmetic  operations  is  dominated  by  0(K(m^  +  m^nQ)  -t-n®) 

Table  1:  Complexity  Estimates  for  Dominant  Steps  in  Method  of  §4 


Equation  Number 

Estimate  of  the  Number 

of  Computation 

of  Arithmetic  Operations 

(21) 

(25) 

(28) 

(30) 

(32) 

(34) 

0(K(m.Q  +  m0no)) 
OlyKmQUQ  +  Uq) 
0(m^n^  +  mg) 
O(m0nl) 

0(m^n^  +  mg) 
O(monl) 

for  all  three  choices  of  P  indicated  in  §2.2. 

To  analyze  the  work  of  computation  (25),  we  let 

svec{BlM~lBkU)  =  svec{J2i  •  u)Bk]) 

=  •U)svec(B^) 

=  £*  £j  4$  (svec(B^  )T svec(U))svec(B^ ) 

=  £i  £j  4j\svec(B[!))svec(Blf))T)svec(U). 

So  matrix  Mo  =  Fq  1  So  +  J2k= l  BkMjf1Bk  in  Mn°Vn°  is  given  by 

K  mo  mo 

Ho  +  4j)(svec{B[!))svec(B^})T),  (36) 

fc=l  i= 1  j= 1 

where  Ho  is  J-ffl£o,  which  depends  on  different  choices  of  symmetrization  for  search 
directions.  The  number  of  arithmetic  operators  in  (36)  is  0(Km^nQ+nfy.  This  completes 
the  proof.  □ 
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If  we  use  a  generic  homogeneous  algorithm  such  as  the  one  in  [14],  then  the  number  of 
arithmetic  operations  required  to  compute  the  search  directions  for  (11,  13)  is  0(mn3  + 
m2n2  +  to3),  where  n  =  (no  +  Km)  and  to  =  (mo  +  Kim).  Setting  m\  =  mo  and 
m  =  no  and  substituting  n  =  (1  +  K)no  and  m  =  (1  +  K)mo)  in  0(mn3  +  m2n 2  +  ?n3), 
we  have  that  the  complexity  of  such  a  generic  method  of  computing  the  search  directions 
is  0(K4(monQ  +  TOqU-q)  +  K3m .g).  This  is  much  larger  than  the  complexity  0{K(m g  + 
mono)  +  no)  obtained  for  the  method  in  §4  when  K  mo  and  K  no- 

It  has  been  shown  in  [11]  that  if  Problem  (11,  13)  has  a  solution,  then  the  KSH 
method  is  globally  convergent.  The  algorithm  finds  an  optimal  solution  or  determines 
that  the  primal-dual  pair  has  no  solution  of  norm  less  than  a  given  number  in  at  most 
0(n1/2L)  iterations,  where  n  is  the  size  of  the  problem  and  L  is  the  logarithm  of  the 
ratio  of  the  initial  error  and  the  tolerance.  So  by  using  the  method  in  §4  for  computing 
the  search  direction  with  KSH  symmetrization,  the  complexity  of  Algorithm  1  in  terms 
of  the  total  number  of  arithmetic  operations  is  0(K'^2 {mfpnj1  -t-m-gn^2)  +  K1/2^).  In 
comparison,  the  short-  and  long-step  decomposition  algorithms  of  Ariyawansa  and  Zhu 
[2]  have  complexities  of  0(K 3/2)  and  0(K2)  respectively  in  terms  of  the  total  number 
of  arithmetic  operations. 

We  note  that  the  efficient  computation  of  the  Schur  computation  matrix  M k  in  (21) 
is  crucial  as  this  is  the  most  expensive  step  in  each  iteration  where  usually  80%  of  the 
total  CPU  time  is  spent  if  the  algorithm  in  [15]  is  used.  However,  in  each  iteration, 
each  Mk  can  be  computed  independently,  and  so  distributed  processing  may  be  used  to 
achieve  substantial  reductions  in  computation  time. 

6  Concluding  Remarks 

In  this  paper,  we  have  given  a  homogeneous  self-dual  algorithm  for  Stochastic  Semidefi- 
nite  Programming.  Being  a  homogeneous  self-dual  algorithm  it  does  not  require  the  user 
to  provide  a  starting  point.  Finding  a  starting  point  may  be  difficult  for  SSDP  prob¬ 
lems.  We  have  also  developed  an  efficient  method  for  calculating  the  search  directions 
which  can  take  advantage  of  parallel  and  distributed  processing.  A  worst-case  bound  on 
the  number  of  arithmetic  operations  for  our  algorithm  has  been  obtained.  This  bound 
shows  that  the  complexity  of  our  algorithm  is  similar  to  that  of  volumetric  decomposition 
algorithms  described  in  [2]. 
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